
library(gdxrrw)
library(ggplot2)
library(dplyr)
library(reshape2)
library(tidyr)
library(grid)
library(RColorBrewer)
library(gridExtra)
library(tableone)
library(cowplot)
library(png)
library(sp)
library(maptools)
library(colorspace)
library(maps)
library(scales)

OrRdPal <- brewer.pal(9, "OrRd")
GnBuPal <- brewer.pal(9, "GnBu")
YlGnBupal <- brewer.pal(9, "YlGnBu")
Redspal <- brewer.pal(9, "Reds")
bluesPal <- brewer.pal(9, "Blues")
pastelpal <- brewer.pal(9, "Pastel1")
dark2pal <- brewer.pal(8, "Dark2")
Set1pal <- brewer.pal(8, "Set1")
alphabet <- c("a","b","c","d","e","f","g","h","i","j","k")

system("gams DataProcess.gms")

MyThemeLine <- theme_bw() +
  theme(
    panel.border=element_rect(fill=NA),
    panel.grid.minor = element_line(color = NA), 
    #    axis.title=element_text(size=5),
    #    axis.text.x = element_text(hjust=1,size = 10, angle = 0),
    axis.line=element_line(colour="black"),
    panel.background=element_rect(fill = "white"),
    #    panel.grid.major=element_line(linetype="dashed",colour="grey",size=0.5),
    panel.grid.major=element_blank(),
    strip.background=element_rect(fill="white", colour="white"),
    strip.text.x = element_text(size=10, colour = "black", angle = 0,face="bold"),
    axis.text.x=element_text(size = 10,angle=45, vjust=0.9, hjust=1, margin = unit(c(t = 0.3, r = 0, b = 0, l = 0), "cm")),
    axis.text.y=element_text(size = 10,margin = unit(c(t = 0, r = 0.5, b = 0, l = 0.5), "cm")),
    legend.text = element_text(size = 10),
    legend.title = element_text(size = 8),
    axis.ticks.length=unit(0.15,"cm")
  )


#-- data load
dir.create("../output/")
dir.create("../output/4paper")
dir.create("../output/gdx")
dir.create("../output/png")
dir.create("../output/txt")

#-file copy from server
ServCopy <- c("AnalysisExpenditure","AnalysisIncome")
for(i in 1:length(ServCopy)){
  v <- c(paste0("//10.249.219.93/sfujimori/Pov/Pov3/output/gdx/",ServCopy[i],".gdx"),paste0("../data/",ServCopy[i],".gdx"))
  file.copy(v[1], v[2], overwrite = TRUE) 
}
ServCopy <- c("Inputdata")
for(i in 1:length(ServCopy)){
  v <- c(paste0("//10.249.219.93/sfujimori/Pov/Pov3/PovertyIncome/data/",ServCopy[i],".gdx"),paste0("../data/",ServCopy[i],".gdx"))
  file.copy(v[1], v[2], overwrite = TRUE) 
}
ServCopy <- c("global_17_iamc","analysis")
for(i in 1:length(ServCopy)){
  v <- c(paste0("//10.249.219.93/sfujimori/Pov/Pov3/PovertyIncome/data/cgeoutput/",ServCopy[i],".gdx"),paste0("../data/",ServCopy[i],".gdx"))
  file.copy(v[1], v[2], overwrite = TRUE) 
}


outputdir <- c("../output/")
linepalette <- c("#FF7F00","#4DAF4A","#377EB8","#E41A1C","#984EA3","#A65628","#F781BF","#8DD3C7","#FB8072","#80B1D3","#FDB462","#B3DE69","#FCCDE5","#D9D9D9","#BC80BD","#CCEBC5","#FFED6F","#7f878f")
landusepalette <- c("#8DD3C7","#FF7F00","#377EB8","#4DAF4A","#A65628")
region <- c("World","ASIA2","CHN","IND","JPN","XSE","XSA")
scenariomap <- read.table("../define/scenariomap.map", sep="\t",header=T, stringsAsFactors=F)
ScenarioMapPov <- read.table("../define/ScenarioMapPov.map", sep="\t",header=T, stringsAsFactors=F)
ScenarioMapDist <- read.table("../define/scenarioDistribtion.txt", sep="\t",header=T, stringsAsFactors=F)

ScenarioMapCGE <- read.table("../define/ScenarioMapCGE.map", sep="\t",header=T, stringsAsFactors=F)
varlist_load <- read.table("../define/all_list.txt", sep="\t",header=F, stringsAsFactors=F) %>% rename(VEMF=V1,VName=V2,Unit=V3)
countrycodelist <- read.table("../define/WDIcountry.set", sep="\t",header=F, stringsAsFactors=F) %>% rename(R=V1,Rname=V2)
tpespalette <- c("Coal|w/o CCS"="#000000","Coal|w/ CCS"="#7f878f","Oil|w/o CCS"="#ff2800","Oil|w/ CCS"="#ffd1d1","Gas|w/o CCS"="#9a0079","Gas|w/ CCS"="#c7b2de","Hydro"="#0041ff","Nuclear"="#663300","Solar"="#b4ebfa","Wind"="#ff9900","Biomass|w/o CCS"="#35a16b","Biomass|w/ CCS"="#cbf266","Geothermal"="#edc58f","Other"="#ffff99",
                 "Solid"=pastelpal[1],"Liquid"=pastelpal[2],"Gas"=pastelpal[3],"Electricity"=pastelpal[4],"Heat"=pastelpal[5],
                 "Build-up"=Set1pal[1],"Cropland (for food)"=Set1pal[2],"Cropland"=Set1pal[2],"Forest"=Set1pal[3],"Pasture"=Set1pal[4],"Energy Crops"=Set1pal[5],"Other Land"=Set1pal[6],"Other Arable Land"=Set1pal[7])
areamap <- read.table("../define/Areafigureorder.txt", sep="\t",header=T, stringsAsFactors=F)
areamappara <- read.table("../define/Area.map", sep="\t",header=T, stringsAsFactors=F)
R5Name <- read.table("../define/R5Name.txt", sep="\t",header=T, stringsAsFactors=F)
ragg <- c("R5OECD90+EU","R5ASIA","R5LAM","R5MAF","R5REF","WLD")
R2Name <- R5Name %>% filter(R %in% c("R5ASIA","R5MAF"))

SceOrder <- read.table("../define/Scenarioorder.txt", sep="\t",header=T, stringsAsFactors=F)

#loading model results
#-- CGE 
CGEload0 <- rgdx.param('../data/global_17_iamc.gdx','IAMC_Template') %>% inner_join(varlist_load) %>% rename("Value"=IAMC_template,R=REMF) %>% 
  left_join(scenariomap,by="SCENARIO") %>% filter(SCENARIO %in% as.vector(scenariomap[,1]))
Giniload <- rgdx.param('../data/Inputdata.gdx','Gini') %>% rename("Value"=Gini,SCENARIO=Ref)%>% mutate(VEMF="Gini") %>% inner_join(varlist_load) %>%  
  left_join(ScenarioMapPov,by="SCENARIO") %>% filter(SCENARIO %in% as.vector(ScenarioMapPov[,1])) %>% select(SCENARIO,R,VEMF,Y,Value,VName,Unit,SceName,SocEco)  
PriceChangeload <- rgdx.param('../data/Inputdata.gdx','PriceChange') %>% rename(Value=PriceChange,SCENARIO=Ref)%>%  left_join(ScenarioMapPov,by="SCENARIO") %>% filter(SCENARIO %in% as.vector(ScenarioMapPov[,1]))   
PovHeadHis <- rgdx.param('../data/Inputdata.gdx','PovertyHeadcount') %>% rename(Value=PovertyHeadcount,SCENARIO=Ref)%>%  left_join(ScenarioMapPov,by="SCENARIO") %>% filter(SCENARIO %in% "WDI")  %>% mutate(Model="Income",Ind="PovHeadC")
PovHeadHis$Y <- as.numeric(levels(PovHeadHis$Y))[PovHeadHis$Y]

CarbTaxRateload <- rgdx.param('../output/gdx/CGEDataprocess.gdx','CTAXGDPRateScenario') %>% rename(Value=CTAXGDPRateScenario,SCENARIO=SCENARIOCGE) %>%  left_join(ScenarioMapCGE,by="SCENARIO") %>% filter(SCENARIO %in% as.vector(ScenarioMapCGE[,1]))   
CarbTaxRateload$Y <- as.numeric(levels(CarbTaxRateload$Y))[CarbTaxRateload$Y]
DifCarbT_PovGapload <- rgdx.param('../output/gdx/CGEDataprocess.gdx','PovGap_CtaxRev_Dif') %>% rename(Value=PovGap_CtaxRev_Dif,SCENARIO=SCENARIOCGE,Ind=".i6") %>%  left_join(ScenarioMapCGE,by="SCENARIO") %>% filter(SCENARIO %in% as.vector(ScenarioMapCGE[,1]))   
DifCarbT_PovGapload$Y <- as.numeric(levels(DifCarbT_PovGapload$Y))[DifCarbT_PovGapload$Y]
CarbTaxload <- rgdx.param('../output/gdx/CGEDataprocess.gdx','CTAXScenario') %>% rename(Value=CTAXScenario,SCENARIO=SCENARIOCGE) %>%  left_join(ScenarioMapCGE,by="SCENARIO") %>% filter(SCENARIO %in% as.vector(ScenarioMapCGE[,1]))   
CarbTaxload$Y <- as.numeric(levels(CarbTaxload$Y))[CarbTaxload$Y]

Giniload$R <- ifelse(Giniload$R=="WLD", "World", as.vector(Giniload$R))
CGEloadfull <- rbind(CGEload0,Giniload)
CGEloadfull$Y <- as.numeric(levels(CGEloadfull$Y))[CGEloadfull$Y]

#-- POverty tool
PovHeadExpenditure <- rgdx.param('../data/AnalysisExpenditure.gdx','PovExp') %>% rename("Value"=PoVExp,SCENARIO=Ref) %>% mutate(Model="Expenditure",Ind="PovHeadC")  
PovHeadIncome <- rgdx.param('../data/AnalysisIncome.gdx','PoV') %>% rename("Value"=PoV,SCENARIO=Ref) %>% mutate(Model="Income",Ind="PovHeadC") 
RoPIncome <- rgdx.param('../data/AnalysisIncome.gdx','RoP') %>% rename("Value"=RoP,SCENARIO=Ref) %>% mutate(Model="Income",Ind="RoP") 
PovGapGDPExpenditure <- rgdx.param('../data/AnalysisExpenditure.gdx','PoVgap_GDPExp') %>% rename("Value"=PoVgap_GDPExp,SCENARIO=Ref) %>% mutate(Model="Expenditure",Ind="PovGapGDP") 
PovGapGDPIncome <- rgdx.param('../data/AnalysisIncome.gdx','Povgap_GDP') %>% rename("Value"=Povgap_GDP,SCENARIO=Ref) %>% mutate(Model="Income",Ind="PovGapGDP") 
PovGapExpenditureAbs <- rgdx.param('../data/AnalysisExpenditure.gdx','PoVgap_absExp') %>% rename("Value"=PoVgap_absExp,SCENARIO=Ref) %>% mutate(Model="Expenditure",Ind="PovGapAbs") 
PovGapIncomeAbs <- rgdx.param('../data/AnalysisIncome.gdx','Povgap_abs') %>% rename("Value"=Povgap_abs,SCENARIO=Ref) %>% mutate(Model="Income",Ind="PovGapAbs") 
PovRateExpenditure <- rgdx.param('../data/AnalysisExpenditure.gdx','RoPExp') %>% rename("Value"=RoPExp,SCENARIO=Ref) %>% mutate(Model="Expenditure",Ind="PovRate")  
PovRateIncome <- rgdx.param('../data/AnalysisIncome.gdx','RoP') %>% rename("Value"=RoP,SCENARIO=Ref) %>% mutate(Model="Income",Ind="PovRate") 
PovGapExpenditure <- rgdx.param('../data/AnalysisExpenditure.gdx','PoVgap_rateExp') %>% rename("Value"=PoVgap_rateExp,SCENARIO=Ref) %>% mutate(Model="Expenditure",Ind="PovGap") 
PovGapIncome <- rgdx.param('../data/AnalysisIncome.gdx','Povgap') %>% rename("Value"=Povgap,SCENARIO=Ref) %>% mutate(Model="Income",Ind="PovGap") 
GDPCap <- rgdx.param('../data/Inputdata.gdx','GDPCap') %>% rename("Value"=GDPCap,SCENARIO=Ref)%>% left_join(ScenarioMapPov,by="SCENARIO") %>% filter(SCENARIO %in% as.vector(ScenarioMapPov[,1]))
GDPCap$Y <- as.numeric(levels(GDPCap$Y))[GDPCap$Y]

PovIndicator <- rbind(PovHeadExpenditure,PovHeadIncome,PovGapGDPExpenditure,PovGapGDPIncome,PovRateExpenditure,PovRateIncome,PovGapIncomeAbs,RoPIncome,PovGapExpenditure,PovGapIncome) %>% left_join(ScenarioMapPov,by="SCENARIO") %>% filter(SCENARIO %in% as.vector(ScenarioMapPov[,1])) %>% rbind(PovHeadHis)
PovIndicator$Y <- as.numeric(levels(PovIndicator$Y))[PovIndicator$Y]
PovIndicatorIncomeBaU <-  filter(PovIndicator,Model=="Income" & SceName=="Baseline") 
PovIndicatorIncomeBaU$Model <- "Expenditure"
PovIndicator0 <- rbind(PovIndicator,PovIndicatorIncomeBaU)

PovIndicatorExemptBase <- PovIndicator0 %>% filter(SocEco=="SSP2" & Model=="Income" & SceName=="Baseline") %>% select(-SocEco) %>% mutate(SocEco="SSP2_exempt")
PovIndicator <- rbind(PovIndicator0,PovIndicatorExemptBase) %>% left_join(SceOrder)

PovHeadC <- filter(PovIndicator,Ind=="PovHeadC")
PovGap <- filter(PovIndicator,Ind=="PovGap")
PovRate <- filter(PovIndicator,Ind=="PovRate")
PovGapAbs <- filter(PovIndicator,Ind=="PovGapAbs")
RoP <- filter(PovIndicator,Ind=="RoP")

#Bar graph for poverty implications of mitigation policy
Indlist <- read.table("../define/IndicatorList.txt", sep="\t",header=F, stringsAsFactors=F) %>% rename(Ind=V1,IndName=V2,Unit=V3,scaler=V4)
plot.mit <- as.list(as.vector(Indlist$Ind))
names(plot.mit) <- as.vector(Indlist$Ind)
for(v in Indlist$Ind){
  ylabel <- Indlist[Indlist$Ind==v,]$IndName
  ylabelunit <- Indlist[Indlist$Ind==v,]$Unit
  YY <- filter(left_join(filter(PovIndicator,Ind %in% Indlist$Ind),Indlist,by="Ind"),Ind==v,Y==2030 & R=="WLD"& TH=="pop_1.9" & SocEco %in% c("SSP2","SSP2_exempt"))
  X2 <- ggplot() +
    geom_bar(data=YY,aes(x=reorder(X=Sceorder,x=SceName),y=Value*scaler,fill=reorder(X=Sceorder,x=SceName)),stat="identity", position="identity")+MyThemeLine +
    xlab("Scenarios") + ylab(paste0(ylabel,ylabelunit)) +scale_fill_manual(values=linepalette)+ theme(legend.title=element_blank()) 
  X3 <-  X2 + facet_wrap(SocEco~ Model,ncol=4) 
  plot.mit[[v]] <- X3
}
ggsave(plot.mit[["PovRate"]], file=paste0("../output/png/Incomeexpendigure_Povertyrate.png"), dpi = 150, width=10, height=4,limitsize=FALSE)
ggsave(plot.mit[["PovHeadC"]], file=paste0("../output/png/Incomeexpendigure_PHC.png"), dpi = 150, width=10, height=4,limitsize=FALSE)

#--- Decomposition of factors
PovIndicatortmp1 <- PovIndicator %>% filter(Model=="Income" & SceName!="Historical" & SocEco=="SSP2") %>% select (-Sceorder,-SCENARIO) %>% spread(key=SceName,value=Value) %>% select(-SocEco)
PovIndicatortmp100 <- PovIndicatortmp1 %>% select(-R,-Y,-TH,-Model,-Ind)
PovIndicatortmp100 <- PovIndicatortmp100 - PovIndicatortmp100$Baseline
PovIndicatortmp1.1 <- gather(cbind(select(PovIndicatortmp1,R,Y,TH,Model,Ind),PovIndicatortmp100),key=SceName,value=IncomeEff,-R,-Y,-TH,-Model,-Ind) %>% left_join(SceOrder)

PovIndicatortmp2 <- PovIndicator %>% filter(Model=="Expenditure" & SceName!="Historical" & SocEco=="SSP2") %>% select (-Sceorder,-SCENARIO) %>% spread(key=SceName,value=Value) %>% select(-SocEco)
PovIndicatortmp200 <- PovIndicatortmp2 %>% select(-R,-Y,-TH,-Model,-Ind)
PovIndicatortmp200 <- PovIndicatortmp200 - PovIndicatortmp200$Baseline
PovIndicatortmp2.1 <- gather(cbind(select(PovIndicatortmp2,R,Y,TH,Model,Ind),PovIndicatortmp200),key=SceName,value=TotalEff,-R,-Y,-TH,-Model,-Ind) %>% left_join(SceOrder) %>% select(-Model)

PovIndicatortmp4 <- PovIndicator %>% filter(Model=="Expenditure")%>% select (-SCENARIO) %>% spread(key=SocEco,value=Value) %>% select(-Model) %>% na.omit()
PovIndicatortmp4$CtaxDirect <- PovIndicatortmp4$SSP2-PovIndicatortmp4$SSP2_exempt
PovIndicatortmp4.1 <- inner_join(inner_join(PovIndicatortmp1.1,PovIndicatortmp4),PovIndicatortmp2.1) 
PovIndicatortmp4.1$PriceEff <- PovIndicatortmp4.1$TotalEff-PovIndicatortmp4.1$IncomeEff-PovIndicatortmp4$CtaxDirect
PovIndicatortmp3 <- PovIndicatortmp4.1 %>% select("R", "Y", "TH", "Ind","SceName","IncomeEff","CtaxDirect","PriceEff") %>% gather(key=Eff,value=value,-R,-Y,-TH,-Ind,-SceName)
FacOrder <- data.frame("Eff"=c("CtaxDirect","IncomeEff","PriceEff"),"EffOrd"=c(3,1,2))

DecomList <- list(1)
for(v in Indlist$Ind){
  for(rr in c("WLD")){
    YY <- filter(left_join(filter(PovIndicatortmp3,Ind %in% Indlist$Ind),Indlist,by="Ind"),Ind==v,Y %in% c(2015,2030,2050) & R==rr & TH=="pop_1.9") %>% left_join(SceOrder,by=c("SceName")) %>% left_join(FacOrder)
    ylabel <- Indlist[Indlist$Ind==v,]$IndName
    ylabelunit <- Indlist[Indlist$Ind==v,]$Unit
    plot1 <- ggplot() + geom_bar(data=filter(YY,value>=0),aes(x=reorder(X=Sceorder,x=SceName),y=value*scaler,fill=reorder(Eff,EffOrd)),stat="identity")+MyThemeLine +
      xlab("Scenarios") + ylab(paste0(ylabel,ylabelunit)) +scale_fill_manual(values=linepalette)+ theme(legend.title=element_blank()) +facet_wrap(Y~ .,ncol=4) 
    ggsave(plot1, file=paste0("../output/png/EffectsAbs",rr,v,".png"), dpi = 150, width=10, height=4,limitsize=FALSE)
    plot3 <-  ggplot() + geom_bar(data=filter(YY,value>=0),aes(x=reorder(X=Sceorder,x=SceName),y=value*scaler,fill=reorder(Eff,EffOrd)),stat="identity",position = "fill")+MyThemeLine +
      xlab("Scenarios") + ylab(paste0(ylabel,ylabelunit)) +scale_fill_manual(values=linepalette)+ theme(legend.title=element_blank()) +facet_wrap(Y~ .,ncol=4) + scale_y_continuous(labels = percent)  
    ggsave(plot3, file=paste0("../output/png/EffectsShare",rr,v,".png"), dpi = 150, width=10, height=4,limitsize=FALSE)
    pp6 <- plot_grid(plot1,plot3,ncol=1) +  
      draw_plot_label(label = alphabet[1:2], size = 12,x = c(0.00,0.00), y = c(1,0.5))
    ggsave(pp6, file=paste0("../output/png/DecompWorld",v,".png"), dpi = 250, width=7, height=7,limitsize=FALSE)

    Decompplot <- list(plot1)
    names(Decompplot) <- c(paste0("DecAbs",v))
    DecomList <- c(DecomList,Decompplot)
    Decompplot <- list(plot3)
    names(Decompplot) <- c(paste0("DecShr",v))
    DecomList <- c(DecomList,Decompplot)
  }
}


#Regional
v <- Indlist$Ind[2]
YY <- filter(left_join(filter(PovIndicatortmp3,Ind %in% Indlist$Ind),Indlist,by="Ind"),Ind==v,Y %in% c(2030) & R %in% R2Name$R & TH=="pop_1.9") %>% left_join(SceOrder,by=c("SceName")) %>% left_join(FacOrder)
ylabel <- Indlist[Indlist$Ind==v,]$IndName
ylabelunit <- Indlist[Indlist$Ind==v,]$Unit
plot1 <- ggplot() + geom_bar(data=filter(YY,value>=0),aes(x=R,y=value*scaler,fill=reorder(Eff,EffOrd)),stat="identity")+MyThemeLine +
  xlab("Scenarios") + ylab(paste0(ylabel,ylabelunit)) +scale_fill_manual(values=linepalette)+ theme(legend.title=element_blank()) +facet_wrap(reorder(X=Sceorder,x=SceName)~ .,ncol=4) 
plot2 <- ggplot() + geom_bar(data=filter(YY,value>=0),aes(x=R,y=value*scaler,fill=reorder(Eff,EffOrd)),stat="identity",position = "fill")+MyThemeLine + scale_y_continuous(labels = percent)+
  xlab("Scenarios") + ylab(paste0(ylabel,ylabelunit)) +scale_fill_manual(values=linepalette)+ theme(legend.title=element_blank()) +facet_wrap(reorder(X=Sceorder,x=SceName)~ .,ncol=4) 

pp6 <- plot_grid(plot1,plot2,ncol=1) +  
  draw_plot_label(label = alphabet[1:2], size = 12,x = c(0.00,0.00), y = c(1,0.5))
ggsave(pp6, file="../output/png/DecompRegion.png", dpi = 250, width=5, height=5,limitsize=FALSE)
Decompplot <- list(plot1)
names(Decompplot) <- c(paste0("DecAbsRegion",v))
DecomList <- c(DecomList,Decompplot)
Decompplot <- list(plot2)
names(Decompplot) <- c(paste0("DecShrRegion",v))
DecomList <- c(DecomList,Decompplot)

for(sce in c("WB2C","1.5C")){
  YY <- filter(left_join(filter(PovIndicatortmp3,Ind %in% Indlist$Ind),Indlist,by="Ind"),Ind==v,Y %in% c(2030) & !(R %in% R5Name$R) & !(R %in% c("WLD","ASIA2")) & TH=="pop_1.9" & SceName==sce) %>% spread(key=Eff,value=value)
  YY$PovHeadC <- YY$CtaxDirect + YY$IncomeEff + YY$PriceEff
  YY <- YY[order(YY$"PovHeadC",decreasing =T),] %>% mutate(Povnum = row_number()) %>% left_join(countrycodelist)
  YY2 <- YY %>% select (-PovHeadC) %>% gather(key=Eff,value=value,"CtaxDirect","IncomeEff","PriceEff")  %>% left_join(FacOrder)
  YY3 <- select(YY,CtaxDirect,IncomeEff,PriceEff,PovHeadC)
  YY3 <- YY3/YY3$PovHeadC
  YY4 <- cbind(select(YY,-CtaxDirect,-IncomeEff,-PriceEff,-PovHeadC),YY3) %>% select (-PovHeadC) %>% gather(key=Eff,value=value,"CtaxDirect","IncomeEff","PriceEff") %>% left_join(FacOrder)
  
  plot1.country <- ggplot() + geom_bar(data=filter(YY2,Povnum<=30),aes(x=reorder(Rname,Povnum),y=value*scaler,fill=reorder(Eff,EffOrd)),stat="identity")+MyThemeLine +
    xlab("countries") + ylab(paste0(ylabel,ylabelunit)) +scale_fill_manual(values=linepalette)+ theme(legend.title=element_blank()) 
  plot2.country <- ggplot() + geom_bar(data=filter(YY4,Povnum<=30),aes(x=reorder(Rname,Povnum),y=value,fill=reorder(Eff,EffOrd)),stat="identity")+MyThemeLine +
    xlab("countries") + ylab(paste0(ylabel,ylabelunit)) +scale_fill_manual(values=linepalette)+ theme(legend.title=element_blank()) 
  pc <- list(plot1.country)
  names(pc) <- c(paste0("abs",sce))
  DecomList <- c(DecomList,pc)
  pc <- list(plot2.country)
  names(pc) <- c(paste0("share",sce))
  DecomList <- c(DecomList,pc)
  
}

#---- Poverty headcount
PovTimeSeriesPlot <- function(v,v2,v3,v4,v5){
  XX <- filter(PovIndicator,Ind==v) %>% filter(SceName %in% v2 & R %in% v3 & TH %in% v4 & Y>=ystart & Y<=yend & Model %in% Modelset & SocEco %in% v5)  %>% 
  rbind(filter(PovIndicator,Ind==v & SCENARIO=="Adj" & R %in% regionset& TH %in% v4 )) %>% select(SCENARIO,R,Y,TH,Value,SceName,Sceorder)
  if(v=="PovHeadC"){XX$Value <- XX$Value/10**6}
#  if(length(unique(XX[XX$TH==varset,]$SCENARIO))==1){linepalette3 <- linepalette[2:4]}else{linepalette3 <- linepalette}
  linepalette3 <- linepalette
  PP <- ggplot() + 
    geom_line(data=XX,aes(x=Y, y = Value , color=reorder(SceName,Sceorder),shape=reorder(SceName,Sceorder)),stat="identity", position="identity") +
    geom_point(data=XX,aes(x=Y, y = Value , color=reorder(SceName,Sceorder),shape=reorder(SceName,Sceorder)),size=3.0,fill="white", position="identity") +
    MyThemeLine+ 
        scale_color_manual(values=c("black",linepalette3))  + 
        scale_shape_manual(values=c(1,-1,5,6,7,8,9,10))+
    xlab("year") + ylab(filter(povtslistmap,povtslist==v)[,"povtslabel"])  +  ggtitle(filter(varlist_load,VEMF==varset)[,"VName"]) +
    annotate("segment",x=-Inf,xend=+Inf,y=0,yend=0,linetype="dashed",color="grey")+ theme(legend.title=element_blank()) 
  outname <- paste0("../output/png/",regionset,"_",v,".png")
  ggsave(PP, file=outname, dpi = 150, width=4, height=4,limitsize=FALSE)
  return(PP)
}
PovRegionalPlot <- function(v,v2,v3,v4,v5){
  XX <- filter(PovIndicator,Ind==v) %>% filter(SceName %in% v2 & R %in% v3 & TH %in% v4 & Y>=ystart & Y<=yend & Model %in% Modelset & SocEco %in% v5)  %>% 
    rbind(filter(PovIndicator,Ind==v & SCENARIO=="Adj" & R %in% regionset& TH %in% v4 )) %>% select(SCENARIO,R,Y,TH,Value,SceName)
  if(v=="PovHeadC"){XX$Value <- XX$Value/10**6}
  #  if(length(unique(XX[XX$TH==varset,]$SCENARIO))==1){linepalette3 <- linepalette[2:4]}else{linepalette3 <- linepalette}
  linepalette3 <- linepalette
  PP <- ggplot() + 
    geom_area(data=filter(XX,SceName!="Historical"),aes(x=Y, y = Value , fill=R),stat="identity") +
    MyThemeLine+ scale_fill_manual(values=linepalette3)  + 
    xlab("year") + ylab(filter(povtslistmap,povtslist==v)[,"povtslabel"])  +  ggtitle(filter(varlist_load,VEMF==varset)[,"VName"]) +
    annotate("segment",x=-Inf,xend=+Inf,y=0,yend=0,linetype="dashed",color="grey")+ theme(legend.title=element_blank())+ facet_wrap(SceName ~ .) 
  PP <- PP+geom_area(data=(filter(XX,SceName=="Historical") %>% select(-SceName)),aes(x=Y, y = Value , fill=R),alpha=0.5,stat="identity") 

  PP1 <- ggplot() + 
    geom_area(data=filter(XX,SceName!="Historical"),aes(x=Y, y = Value , fill=R),stat="identity",position = "fill") +
    geom_area(data=(filter(XX,SceName=="Historical") %>% select(-SceName)),aes(x=Y, y = Value , fill=R),alpha=0.5,stat="identity",position = "fill") +
    MyThemeLine+ scale_fill_manual(values=linepalette3)  + 
    xlab("year") + ylab(filter(povtslistmap,povtslist==v)[,"povtslabel"])  +  ggtitle(filter(varlist_load,VEMF==varset)[,"VName"]) +
    annotate("segment",x=-Inf,xend=+Inf,y=0,yend=0,linetype="dashed",color="grey")+ theme(legend.title=element_blank())+ facet_wrap(SceName ~ .)
  
  PPOut <- as.list(c("abs","share"))
  names(PPOut) <- as.vector(c("abs","share"))
  PPOut[["abs"]] <- PP
  PPOut[["share"]] <- PP1
  return(PPOut)
}
scenarioset <- c("Baseline","WB2C","1.5C","NDC")
Soclist<- c("SSP2")
regionset <- c("WLD")
varset <- c("pop_1.9")
yend <- 2050
ystart <- 2010
Modelset <- c("Income")
povtslist <- c("PovHeadC","PovGapGDP","PovGap","PovRate","PovGapCountry","PovGapTaxCountry","PovGapAbsTax","PovGapTaxGDPWorld")
povtslabel <- c("Poverty headcount (million)","Poverty gap per GDP (-)","Poverty gap (-)","Poverty rate (-)","PovGapCountry","PovGapTaxCountry","PovGapAbsTax","PovGapTaxGDPWorld")
povtslistmap <- data.frame(povtslist,povtslabel)
plot.povts <- as.list(povtslist)
names(plot.povts) <- as.vector(povtslist)

Modelset <- c("Expenditure")
for(vv in povtslist[1:4]){
  plot.povts[[vv]] <- PovTimeSeriesPlot(vv,v2=scenarioset,v3=regionset,v4=varset,v5=Soclist)
}
#--- regional figure
regionset <- c("R5ASIA","R5OECD90+EU","R5MAF","R5REF","R5LAM")
scenarioset <- c("Baseline","WB2C")
for(vv in povtslist[1]){
  plot.reg <- PovRegionalPlot(vv,v2=scenarioset,v3=regionset,v4=varset,v5=Soclist)
}

#----price figures
#Settting parameters for price information
scenarioset <- c("1.5C", "WB2C", "NDC")
regionset <- c("WLD")
ylist <- c("2030","2050")
pricelist <- c("Relative","2030","2030Exempt")
plot.price <- as.list(pricelist) 
names(plot.price) <- pricelist
PricePlot <- function(XX){
  linepalette3 <- linepalette
  XX$SceName = factor(XX$SceName, levels = c( "NDC", "WB2C", "1.5C"))
  plotreturn <- ggplot() + 
    geom_bar(data=XX,aes(y=100*(Value-1), x = I , fill=SceName),stat="identity", position="identity")+ coord_flip() +
    MyThemeLine+ scale_fill_manual(values=linepalette[2:5])  +
    xlab("Goods") + ylab("Price index percentage changes(%)")  +  ggtitle("") + theme(legend.position="none") +
    annotate("segment",y=0,yend=0,x=-Inf,xend=+Inf,linetype="dashed",color="grey")+ theme(legend.title=element_blank()) 
  plotreturn <- plotreturn + facet_wrap(SceName~Y,ncol=2)
  return(plotreturn)
}
PricePlot2 <- function(XX){
  linepalette3 <- linepalette
  XX$SceName = factor(XX$SceName, levels = c("Baseline", "NDC", "WB2C", "1.5C"))
  plotreturn <- ggplot() + 
    geom_bar(data=XX[XX$Y==2030,],aes(y=100*(Value-1), x = SceName , fill=SceName),stat="identity", position="identity")+ coord_flip() +
    MyThemeLine+ scale_fill_manual(values=linepalette[2:5])  +
    xlab("Goods") + ylab("Price index percentage changes(%)")  +  ggtitle("") + theme(legend.position="bottom") +
    annotate("segment",y=0,yend=0,x=-Inf,xend=+Inf,linetype="dashed",color="grey")+ theme(legend.title=element_blank()) 
  plotreturn <- plotreturn + facet_wrap(~I,ncol=3)
  return(plotreturn)
}
PriceChangetmp0 <- PriceChangeload %>% filter(SocEco=="SSP2" & SceName=="Baseline")
PriceChangetmp0$SocEco ="SSP2_exempt"
PriceChangetmp <- rbind(PriceChangeload,PriceChangetmp0) %>% select(-SCENARIO) %>% spread(key=SceName,value=Value)
PriceChangetmp1 <- PriceChangetmp %>% select(-R,-Y,-I,-SocEco) 
PriceChangetmp1 <- PriceChangetmp1 / PriceChangetmp1$Baseline
PriceChangedata1 <- gather(cbind(select(PriceChangetmp,R,Y,I,SocEco),PriceChangetmp1),key=SceName,value=Value,-R,-Y,-I,-SocEco) %>% left_join(SceOrder)

data4plot <- PriceChangedata1 %>% filter(SceName %in% scenarioset & R %in% regionset & Y %in% ylist & SocEco %in% c("SSP2"))
plot.price[[1]] <- PricePlot(data4plot)
ggsave(plot.price[[1]], file=paste0("../output/png/",regionset,"_PriceChange.png"), dpi = 150, width=7, height=9,limitsize=FALSE)
plot.price[[2]] <- PricePlot2(data4plot)
ggsave(plot.price[[2]], file=paste0("../output/png/",regionset,"_PriceChange2030.png"), dpi = 150, width=12, height=9,limitsize=FALSE)
data4plot <- PriceChangedata1 %>% filter(SceName %in% scenarioset & R %in% regionset & Y %in% ylist & SocEco %in% c("SSP2","SSP2_exempt"))
plot.price[[3]] <- PricePlot2(data4plot)
plot.price[[3]] <- plot.price[[3]] + facet_wrap(SocEco~I,ncol=3)
outname <- paste0("../output/png/",regionset,"_PriceChange2030Exempt.png")
ggsave(plot.price[[3]], file=outname, dpi = 150, width=12, height=9,limitsize=FALSE)



#---- Distribution
Distribution_income_load    <- rgdx.param('../data/AnalysisIncome.gdx','ISDist10') %>% rename(Value="ISDist10",IS=IS10,Y=Y10) %>% mutate(Model="Income",Step="10")
Distribution_expenditure_load    <- rgdx.param('../data/AnalysisExpenditure.gdx','ISDist10_exp') %>% rename(Value="ISDist10_exp",Y=Y10,IS=IS10) %>% mutate(Model="Expenditure",Step="10")
Distribution_load2 <- rbind(Distribution_income_load,Distribution_expenditure_load)
Distribution_load2$IS <- as.numeric(levels(Distribution_load2$IS))[Distribution_load2$IS]
Distribution_load2$Y <- as.numeric(levels(Distribution_load2$Y))[Distribution_load2$Y]
Distribution_load3 <-Distribution_load2 %>% left_join(ScenarioMapDist,by="Ref")

#MuExp <- rgdx.param('../data/AnalysisExpenditure.gdx','Mu_exp') %>% mutate(Model="Expenditure",Var="Mu") %>% rename("Value"=Mu_exp,SCENARIO=Ref)%>% left_join(ScenarioMapPov,by="SCENARIO")
#SigmaExp <- rgdx.param('../data/AnalysisExpenditure.gdx','sigma_exp') %>% mutate(Model="Expenditure",Var="Sigma") %>% rename("Value"=sigma_exp,SCENARIO=Ref)%>% left_join(ScenarioMapPov,by="SCENARIO")
#x <-  seq(10, 10000, length=1000)
#Mulog <- filter(rbind(SigmaExp,MuExp),Y==2030 & SceName=="WB2C" & SocEco=="SSP2") %>% spread(key=Var,value=Value)
#fx <- dlnorm(x=x, meanlog = Mu, sdlog=Sigma)
#data <-data.frame(x=x,y=fx)

Ref_list2 <- c("SSP2_Baseline","SSP2_2deg","SSP2_1p5deg")
Ylist3050 <- c(2030,2050)
distlist <- ragg
plot.dist <- as.list(distlist)
names(plot.dist) <- ragg
plot.dist.wrl <- as.list(Ylist3050)
names(plot.dist.wrl) <- Ylist3050
xupper <-5000
for(r in ragg){
  for(ie in c("Expenditure")){
    #Reference across years
    distplotdata <- Distribution_load2 %>% filter(Model==ie & Ref=="SSP2_Baseline" & Step=="10" & Y<=2050 & R==r )
    plot.dist[[r]] <- ggplot()+
      geom_line(data=distplotdata[distplotdata$IS<=xupper/100,],aes(x=IS*100, y=Value/1000,color=interaction(Y)),stat="identity") +
      MyThemeLine + theme(legend.title=element_blank())+
      xlab("Income or Expenditure per capita (2011International$PPP/cap)") + ylab("Population (Thousand)")  +
      annotate("segment",x=365*1.9,xend=365*1.9,y=0,yend=+Inf,linetype="solid",color="grey") 
    ggsave(plot.dist[[r]], file=paste0("../output/png/Distribution",ie,r,"_Year.png"), dpi = 250, width=5, height=4,limitsize=FALSE)
  }
}

#For 2030 across scenarios
for(r in c("WLD")){
  for(ie in c("Expenditure")){
    for(y in c(2030,2050)){
      distplotdata <- Distribution_load3 %>% filter(Model==ie & Ref %in% Ref_list2 & Step=="10" & Y==y & R==r )
      plot.dist.wrl[[y]] <- ggplot()+
        geom_line(data=distplotdata[distplotdata$IS<=xupper/100,],aes(x=IS*100, y=Value/1000,color=SceName),stat="identity") +
        MyThemeLine + theme(legend.title=element_blank())+
        xlab("Income or Expenditure per capita (2011International$PPP/cap)") + ylab("Population (Thousand)")  +
        annotate("segment",x=365*1.9,xend=365*1.9,y=0,yend=+Inf,linetype="solid",color="grey") 
      ggsave(plot.dist.wrl[[y]], file=paste0("../output/png/Distribution",ie,r,y,"_scenarios.png"), dpi = 250, width=5, height=4,limitsize=FALSE)
    }
  }
}


#---- World map
world <- map_data("world")
WorldMapRegion <- unique(world$region)
write.table(WorldMapRegion, file = "../define/WorldMapRegion.txt", append = FALSE, row.names=FALSE, quote = TRUE, sep = ",")
WorldMapRegionMap <- read.table("../define/WorldMapRegion.Map", sep="\t",header=T)

worldmapfunction <- function(X1,X3){
  SDNData <- X1[X1$R=="SDN",]
  SDNData$region <- "South Sudan"
  X2 <- rbind(X1,SDNData)
  WorldMapdata <- left_join(world,X2,by="region")
  g1 <- ggplot(WorldMapdata,aes(x = long, y = lat, group=group, fill = Value)) + geom_polygon(size = 0.25) +
    geom_path(aes(x = long, y = lat, group = group),colour = "gray40", size = 0.01) +xlab("")+ylab("")   +
    scale_fill_gradient2(low = X3[1], mid = X3[2], high = X3[3], na.value="white", midpoint = as.numeric(X3[4]), name=X3[5]
                         #  breaks=c(maxv/6,maxv*3/6,maxv*5/6),labels=c("Low","Medium","High"),
    ) + coord_equal() + scale_x_continuous(breaks=NULL)   +MyThemeLine+
    theme(panel.border = element_rect(colour = "white", fill=NA, size=0),axis.text.x = element_blank(),
          axis.line = element_line(colour="White"),axis.text.y = element_blank(),axis.ticks = element_blank())+ ylim(-65,80)+ xlim(-160,160) + 
    theme(plot.margin = margin(0, 0, 0, 0, "cm"),
          legend.key.width = unit(1, 'cm'), #change legend key size
#          legend.key.height = unit(0.1, 'cm'), #change legend key height
          legend.key.size = unit(1, 'cm'), #change legend key width
          legend.title = element_text(size=14), #change legend title font size
          legend.text = element_text(size=14)) #change legend text font size
  return(g1)  
}

scenariolist <- c("1.5C","NDC","WB2C","Baseline")
yearlist <- c(2010,2030)
jj <- 1
Indlistpov <- c("PovHeadC","PovRate")
Indlistleg <- c("Population at extreme poverty\n(million)","Population at extreme poverty rate (-)")
Indlistmap <- data.frame(Indlistpov,Indlistleg)
WMlist <- list(1)
for(i in Indlistmap$Indlistpov){
  for(yr in yearlist){
    for(v in scenariolist){
      jj = jj +1
      YY <- inner_join(filter(filter(PovIndicator,Ind==i),TH=="pop_1.9" & Y==yr & SceName==v & Model=="Income" & SocEco=="SSP2"),WorldMapRegionMap,by="R")
      if(i=="PovHeadC"){YY$Value <- YY$Value/10**6}
      mappara <- c("white",OrRdPal[4],OrRdPal[7],max(YY$Value)/2,as.vector(filter(Indlistmap,Indlistpov==i)[,"Indlistleg"]))
      WMPovHead <- list(worldmapfunction(YY,mappara) +theme(legend.position="bottom"))
      names(WMPovHead) <- c(paste0(i,v,yr))
      WMlist <- c(WMlist,WMPovHead)
      outname <- paste0("../output/png/Worldmap",i,v,yr,".png")
      ggsave(WMlist[[jj]], file=outname, dpi = 250, width=10, height=6,limitsize=FALSE)
      #Take difference
    }
    YY1 <- spread(inner_join(filter(filter(PovIndicator,Ind==i),TH=="pop_1.9" & Y==yr & Model=="Expenditure" & SocEco=="SSP2"),WorldMapRegionMap,by="R") %>% select(-SCENARIO,-Sceorder),key=SceName,value=Value) 
    YY1$'1.5C' <- YY1$'1.5C'-YY1$'Baseline'
    YY1$'WB2C' <- YY1$'WB2C'-YY1$'Baseline'
    YY1$'NDC' <- YY1$'NDC'-YY1$'Baseline'
    YY2 <- select(YY1,-Baseline) %>% gather(key=SceName,value=Value,'1.5C',WB2C,NDC) %>% left_join(SceOrder)
    for(v in scenariolist){
      if(v!="Baseline" & yr!=2010){
        jj = jj +1
        YY <- filter(YY2, SceName==v)
        if(i=="PovHeadC"){YY$Value <- YY$Value/10**6}
        mappara <- c("white",bluesPal[4],bluesPal[7],max(YY$Value)/2,as.vector(filter(Indlistmap,Indlistpov==i)[,"Indlistleg"]))
        WMPovHead <- list(worldmapfunction(YY,mappara) +theme(legend.position="bottom"))
        names(WMPovHead) <- c(paste0("Dif",i,v,yr))
        WMlist <- c(WMlist,WMPovHead)
        outname <- paste0("../output/png/Worldmap_Dif_",i,v,yr,".png")
        ggsave(WMlist[[jj]], file=outname, dpi = 250, width=10, height=6,limitsize=FALSE)
      }
    }
  }
}

#---- CO2 emissions and poverty comparison
EmissionLoad <- rgdx.param('../data/Emission.gdx','EMI') %>% filter(g=="CO2exLU") %>% rename(R=i_all,Y=t)
EmiPov0 <- inner_join((filter(PovIndicator,Y=="2030" & SceName=="WB2C" & Model=="Expenditure" & SocEco=="SSP2" & TH=="pop_1.9" & Ind %in% c("PovHeadC","PovRate")) %>% spread(key=Ind,value=Value) %>% select(R,PovHeadC,PovRate)),
          (filter(EmissionLoad,Y=="2014") %>% select(R,EMI))) %>% filter(!(R %in% ragg))
EmiPov0.1 <- EmiPov0[order(EmiPov0$"PovRate",decreasing =T),] %>% mutate(Povnum = row_number()) %>% left_join(countrycodelist)
totCO2 <- sum(as.vector(EmiPov0.1$EMI))
totPov <- sum(as.vector(EmiPov0.1$PovHeadC))
EmiPov0.1 <- within(EmiPov0.1, cumCO2 <- cumsum(EMI))
EmiPov0.1 <- within(EmiPov0.1, cumpov <- cumsum(PovHeadC))
EmiPov0.1$cumCO2pct <- EmiPov0.1$cumCO2/totCO2
EmiPov0.1$cumPovpct <- EmiPov0.1$cumpov/totPov
EmiPov0.2 <- EmiPov0.1 %>% select (-PovHeadC,-EMI,-cumCO2,-cumpov) %>% rename("Cumulative CO2 emissions"=cumCO2pct,"Cumulative poverty headcount"=cumPovpct,"Poverty rate"=PovRate) %>% gather(key=Indicator,value=value,"Cumulative CO2 emissions","Cumulative poverty headcount","Poverty rate") 
CO2_PovComp <- ggplot() +
  geom_point(data=filter(EmiPov0.2, Povnum<=50),aes(x=reorder(Rname,Povnum),y = value*100,color=Indicator),stat="identity", position="identity") +
  geom_line(data=filter(EmiPov0.2, Povnum<=50),aes(x=reorder(Rname,Povnum),y = value*100,color=Indicator,group=Indicator),stat="identity", position="identity") +
  MyThemeLine + 
  scale_y_continuous(sec.axis = sec_axis(~ (.) , name = expression("Cumulative percentage of CO2 emissions \n or poverty headcount[%]")),limits = c(0, 100)) +
  xlab("Country") +ylab("Poverty rate (%)") + theme(legend.position="bottom") + theme(plot.margin = grid::unit(c(1, 3, 0.5, 1), "lines"))

ggsave(CO2_PovComp, file=paste0("../output/png/CO2Country.png"), dpi = 250, width=10, height=6,limitsize=FALSE)

#---CGE output
#Setting parameters for selecting CGE output conditions
nalist <- c("POP","GDP_per_cap","GDP_PPP","Gini","Emi_CO2","Prc_Car","Pol_Cos_GDP_Los_rat","Emi_Kyo_Gas","Frc","Tem_Glo_Mea","Prm_Ene","Fin_Ene")
plot.cge <- as.list(nalist)
names(plot.cge) <- nalist
plot.cge2 <- as.list(nalist)
names(plot.cge2) <- nalist
plot.cge3 <- as.list(nalist)
names(plot.cge3) <- nalist
scenarioset <- c("1.5C", "WB2C","NDC","Baseline")
regionset <- c("World","WLD")
regionaggset <- c("R5ASIA","R5OECD90+EU","R5MAF","R5REF","R5LAM")
yend <- 2050
ystart <- 2010
FigList <- list(1)

for(i in nalist){
  for(r in c("World",regionaggset,"IND")){
    varset <- i
    data4plot <- CGEloadfull %>% filter(SceName %in% scenarioset & R %in% r & VEMF %in% varset & Y<= yend  & Y>=ystart & SocEco=="SSP2")%>% left_join(SceOrder)
    if(nrow(data4plot)>=1){
      if(length(unique(data4plot[data4plot$VEMF==varset,]$SCENARIO))==3){linepalette3 <- linepalette[2:4]}else{linepalette3 <- linepalette}
      ylabel <- filter(varlist_load,VEMF==varset)
      plot.cge <- ggplot() + 
        geom_point(data=data4plot,aes(x=Y, y = Value , color=reorder(X=Sceorder,x=SceName)),shape=1,size=3.0,fill="white", position="identity") +
        geom_line(data=data4plot,aes(x=Y, y = Value , color=reorder(X=Sceorder,x=SceName)),stat="identity", position="identity")  +
        MyThemeLine+ scale_color_manual(values=linepalette3)  +
        xlab("Year") + ylab(paste0(ylabel$VName," (",ylabel$Unit,")"))  +  ggtitle(filter(varlist_load,VEMF==varset)[,"VName"]) +
        annotate("segment",x=ystart,xend=yend,y=0,yend=0,linetype="dashed",color="grey")+ theme(legend.title=element_blank()) 
      if(r=="World"){ggsave(plot.cge, file=paste0("../output/png/World_",varset,".png"), dpi = 150, width=4, height=4,limitsize=FALSE)}
      plot.cgeList <- list(plot.cge)
      names(plot.cgeList) <- paste0(r,i)
      FigList <- c(FigList,plot.cgeList)
    }
  }
  
  data4plot <- CGEloadfull %>% filter(SceName %in% c("Baseline") & R %in% regionaggset & VEMF %in% varset & Y>=ystart & Y<=yend & SocEco=="SSP2") %>% inner_join(R5Name,by="R")
  plot.cge2[[i]] <- ggplot() + 
    geom_line(data=data4plot,aes(x=Y, y = Value , color=Rname),stat="identity", position="identity") +
    geom_point(data=data4plot,aes(x=Y, y = Value , color=Rname),shape=1,size=3.0,fill="white", position="identity") +
    MyThemeLine+ scale_color_manual(values=linepalette3)  + theme(legend.position = "none") +
    xlab("Year") + ylab(paste0(filter(varlist_load,VEMF==varset)[,"VName"]," (",filter(varlist_load,VEMF==varset)[,"Unit"],")"))  +  ggtitle(filter(varlist_load,VEMF==varset)[,"VName"]) +
    annotate("segment",x=ystart,xend=yend,y=0,yend=0,linetype="dashed",color="grey")+ theme(legend.title=element_blank()) 
  outname <- paste0("../output/png/region_",varset,".png")
  ggsave(plot.cge2[[i]], file=outname, dpi = 150, width=4, height=4,limitsize=FALSE)

  data4plot <- CGEloadfull %>% filter(SceName %in% scenarioset & R=="World" & VEMF %in% varset & Y>=ystart & Y<=2100 & SocEco=="SSP2")%>% left_join(SceOrder)
  plot.cge3[[i]] <- ggplot() + 
    geom_line(data=data4plot,aes(x=Y, y = Value , color=reorder(X=Sceorder,x=SceName)),stat="identity", position="identity") +
    geom_point(data=data4plot,aes(x=Y, y = Value , color=reorder(X=Sceorder,x=SceName)),shape=1,size=3.0,fill="white", position="identity") +
    MyThemeLine+ scale_color_manual(values=linepalette3)   +
    xlab("Year") + ylab(paste0(filter(varlist_load,VEMF==varset)[,"VName"]," (",filter(varlist_load,VEMF==varset)[,"Unit"],")"))  +  ggtitle(filter(varlist_load,VEMF==varset)[,"VName"]) +
    annotate("segment",x=ystart,xend=2100,y=0,yend=0,linetype="dashed",color="grey")+ theme(legend.title=element_blank()) 
  outname <- paste0("../output/png/",varset,"_Scenario.png")
  ggsave(plot.cge3[[i]], file=outname, dpi = 150, width=4, height=4,limitsize=FALSE)
  
}
## specific figure
data4plot <- CGEloadfull %>% filter(SceName %in% scenarioset & !(R %in% R5Name$R) & VEMF=="Pol_Cos_GDP_Los_rat" & Y==2030 & SocEco=="SSP2")%>% left_join(SceOrder)
ylabel <- filter(varlist_load,VEMF=="Pol_Cos_GDP_Los_rat")
plot.cge.spe1 <- ggplot() + 
  geom_point(data=data4plot,aes(x=R, y = Value , color=reorder(X=Sceorder,x=SceName)),shape=1,size=3.0,fill="white", position="identity") +
  MyThemeLine+ scale_color_manual(values=linepalette[2:4])  +
  xlab("Country") + ylab(paste0(ylabel$VName," (",ylabel$Unit,")"))  +  ggtitle(filter(varlist_load,VEMF==varset)[,"VName"]) 


##merge plots
source("mergeplot.r")
